Layered transition metal dichalcogenides (TMDCs) have attracted extensive interests and been intensively studied in recent years.^[1] Many monolayer TMDCs have been chemical synthesized or mechanically exfoliated successfully from their bulk phases. Crystal structures and electronic properties of TMDCs can be tuned in a wide range by doping, defect structure, adsorption, strain, or stress. Monolayer and few-layered TMDCs show promising properties and performances, and are widely used in catalysis, sensing, optoelectronics, nano-electronic and optical devices.^[2–8]

In addition to widespread interests in the TMDCs for their promising applications, niobium dichalcogenides NbX₂ (X = S, Se, Te) have been widely studied for a long time because of their superconducting behavior and charge density wave (CDW) order.^[9–11] The 2H-NbS₂ is isoelectronic and isostructural with 2H-NbSe₂, showing similar superconducting transition temperature. The CDW instability in NbSe₂ and its possible interaction with superconductivity have been intensive investigated. However, the research progress of CDW in 2H-NbS₂ is stagnant, which has been attributed to the anharmonic suppression.^[12] The nature of the superconducting gap and the pairing mechanism in 2H-NbS₂ are elucidated by ab initio anisotropic Migdal–Eliashberg theory, and the system is proposed to be on the verge of a latent CDW instability.^[13] With the advent of two-dimensional (2D) materials, synthesis of ultrathin, large-size monolayer and few-layered NbS₂ becomes more and more mature.^[14–17] Recently, an enhanced tendency toward CDW order in the 2D limit is observed in NbS₂.^[18–21]

The crystal structure, topological state, and transport properties of the layered materials can be tuned by strain, stress, and pressure, resulting in interesting phenomena and properties such as topological transition,^[22] metal–insulator transition,^[23–25] and ferroelastic phase transition.^[26] As shown in Fig. 1, three polytypes of NbS₂ have been successfully synthesized and feature with one-layer trigonal (1T), two-layer hexagonal (2H), and three-layer rhombohedral (3R) structures.^[27–31] High-pressure experimental investigations indicate that the structure of 3R-NbS₂ is highly anisotropic under compression,^[32,33] and the superconducting transition temperature and upper critical field of 2H-NbS₂ increase simultaneously along with the increasing pressure.^[34] Density functional theory (DFT) calculations predict that a tensile strain can switch the magnetic states in monolayer NbSe₂ and NbS₂.^[35,36] In addition, the structural phase transition and superconducting properties of NbS₂ under high pressure have been studied recently by first-principles calculations, where the 2H phase transforms to a new tetragonal I4/mmm structure at 26 GPa.^[37] The I4/mmm phase [Fig. 1(d)] is predicted to be mechanically and dynamically stable at high pressure by calculating the elastic constants and phonon dispersion, and exhibits a higher superconducting temperature than the 2H phase. Singh et al. discovered an increasing number of conduction bands crossing the Fermi level and higher conductivity in 2H-NbS₂ by applying uniaxial strain along the c-axis.^[38] As shown in Fig. 1, these experimentally already-synthesized phases (1T, 2H, and 3R structures) and theoretically predicted I4/mmm phase crystallize in layered structure. Similar to other layered TMDCs materials, the weak interlayer van der Waals (vdW) interactions should play a crucial role in bonding together the 2D layers in the bulk. Obviously, it is unreasonable to ignore the interlayer vdW interactions in previous theoretical studies on NbS₂ to explore the structural transitions and electronic structures under high pressure or strain.^[37,38]

Fig. 1.

	Figure Option ViewDownloadNew Window
	Fig. 1. The crystal structures of NbS2 with space group (a) (1T), (b) P63/mmc (2H), (c) R3m (3R), and (d) I4/mmm. Big (red) and small (blue) balls represent Nb and S atoms, respectively. The solid line denotes the unit cell. The polyhedra highlight the coordination environments of the Nb ions.

In the present study, we want to revisit the crystal structure and electronic structure of NbS₂ and probe the structural stability of bulk NbS₂ under high pressure by first-principles simulation including vdW corrections. Theoretically calculated results are consistent with the available experimental results of the 2H and 3R phases of NbS₂, indicating that the vdW corrections indeed play a crucial role in the layered crystallographic structures, and improve the description of the electronic states around the Fermi level. Furthermore, we discover that the lattice parameters of the I4/mmm phase remarkably change at 32 GPa, accompanied with abrupt shortening of the interlayer distances and increasing coordination number. The structural transition from layered structure to three-dimensional (3D) bulk in I4/mmm phase can be achieved by imposing a tensile strain in ab-plane, which is originated from the interlayer bonding effects. Our reported results are expected to attract extensive interests from both of the theoretical and experimental communities in the hectic field of TMDCs.

We performed the calculations within the framework of DFT with the projector augmented wave (PAW) method^[39,40] implemented in the Vienna ab initio simulation package (VASP) code.^[41] Generalized gradient approximation (GGA) with the Perdew–Burke–Ernzerhof (PBE) parameterization was used for the exchange–correlation potential.^[42] For geometry optimization and electronic properties calculations, the optimized exchange vdW functional (opt-B88-vdW)^[43] was applied to take into account the vdW interactions. The kinetic-energy cutoff was set up to 520 eV. In order to ensure the accuracy of the structural optimization and allow the atoms to move smoothly, the symmetries were canceled during the relaxing process. Owing to the presence of heavy 4d-electrons in Nb, spin–orbit coupling (SOC) was turned on when calculating the electronic properties. We first confirmed the credibility of the calculations by optimizing the crystal structures of the experimentally already-synthesized phases (1T, 2H, and 3R structures) at ambient condition (0 GPa), then gradually increased the hydrostatic pressure from 0 GPa to 50 GPa with an interval of 5 GPa to simulate the structural responses under compression. The k-points distance was fixed to 0.2 Å⁻¹ in the reciprocal space. The accuracies of the target pressure and energy convergence for all structure optimizations were better than 0.005 GPa and 10⁻⁶ eV, respectively.

In order to examine the influences of the vdW interactions, the crystal structures of the experimentally already-synthesized phases (1T, 2H, and 3R structures) are first optimized at ambient condition (0 GPa) with several vdW-corrected functionals.

As shown in Table 1, the lattice constants obtained by the opt-B88-vdW correction method^[43] agree better with the experimental data^[12,27–30] than the previous reported theoretical results calculated by local density approximate (LDA) method.^[37] We note that although employed the same GGA-PBE functional, even PBE hybrid functional (PBE0),^[44] the calculated unit cell parameters obtained in the present work are closer to the experiment data, implying that the vdW corrections are indeed very important to describe the layered structure of the NbS₂ material. Therefore, the vdW interactions undoubtedly should be taken into account to simulate the structural transitions under high pressure.

Table 1.

Table 1.

Table 1. Comparison of experimental data with theoretically calculated lattice constants a, c and unit cell volume V for the experimentally already-synthesized NbS2 (1T, 2H, and 3R structures) at ambient condition. . Structure Reference a/Å c/Å c/a V/Å3 1T Exp.[27] 3.420 5.938 1.736 60.16 Theo.[37] 3.253 5.341 1.642 48.95 This work 3.384 5.944 1.757 58.96 2H Exp.[12] 3.330 11.950 3.589 114.76 Exp.[27] 3.418 11.860 3.470 119.99 Exp.[28] 3.310 11.890 3.592 112.82 Theo.[37] 3.287 11.421 3.475 106.86 Theo.[44] 3.332 12.106 3.633 116.39 Theo.[44] 3.313 12.074 3.644 114.76 This work 3.358 12.049 3.588 117.66 3R Exp.[28] 3.330 17.810 5.348 171.03 Exp.[29] 3.335 17.834 5.348 171.78 Exp.[30] 3.330 17.918 5.381 172.07 Theo.[37] 3.286 17.577 5.349 164.37 This work 3.357 18.357 5.468 179.10

Table 1. Comparison of experimental data with theoretically calculated lattice constants a, c and unit cell volume V for the experimentally already-synthesized NbS2 (1T, 2H, and 3R structures) at ambient condition. .

Based on the optimized structures, the GGA with opt-B88-vdW correction method is further verified by calculating the band structures of the two common polytypes of NbS₂ (2H and 3R phases). As shown in Fig. 2, the electronic states around the Fermi level for 2H-NbS₂ phase are in good agreement with the electronic structure calculated within GW many-body interaction corrections, where the S p_z and Nb d_z² states are pushed apart at the Γ point.^[45] In contrast to the three bands crossing the Fermi level calculated by LDA or PBE,^[13,44,45] there are only two bands crossing the Fermi level at Γ and K points [Fig. 2(a)], emerging from the bilayer splitting of the Nd d bands. The band structures crossing the Fermi level coincide with the experimental energy and momentum distribution curves.^[46] Especially, the GW calculation results need a rigid shift of the Fermi level to bring in excellent agreement with the experiments.^[45] The present calculations find that the hole pockets centered at Γ (K) predominantly arise from Nb d_z² (d_xy/d_{x² − y²}), and the band dispersions along the Γ–K direction are consistent well with the angle-resolved photoemission spectroscopy (ARPES) experimental results of 2H-NbS₂.^[46] In addition, the obtained band structure of 3R-NbS₂ successfully reproduces the previous theoretical results using various exchange–correlation functionals.^[47,48] The electronic states around Fermi level show similar dispersion characteristics and compositions to those of the 2H phase, arising from two and three repeat trigonal prismatic layers for 2H and 3R NbS₂ [Figs. 1(b) and 1(c)], respectively. SOC only has minor impact on the electronic structures, except that it leads to band splitting around the K point.^[45,48] The present theoretically calculated results indicate that the vdW corrections really play a critical role in bonding together the layered structures. Furthermore, the vdW functional improves the description of the electronic structure through their effects on the crystal structure, giving better agreement with the available experimental results of the 2H and 3R phases of NbS₂.^[49]

Fig. 2.

	Figure Option ViewDownloadNew Window
	Fig. 2. Calculated band structures (solid black lines) and the projection of the electronic states on the level of GGA with opt-B88-vdW correction for the two common polytypes of NbS2: (a) 2H and (b) 3R phases. The orbital contributions are proportional to the sizes of the colored balls, as indicated in the legend: S pz (blue), Nb dz2 (green), and Nb dxy/dx2 − y2 (red).

Figure 3 presents the volume evolutions with compression, which are fitted by the Birch–Murnaghan isothermal equation of state (BM-EOS).^[50] The fitted equilibrium volumes (V₀) are consistent well with the optimized results and experimental data at ambient condition (Table 1). Attributed to the similar layered structure of the experimentally already-synthesized structures (1T, 2H, and 3R phases), their bulk moduli (B₀), the pressure derivative of the bulk modulus (), and volumes at zero pressure (V₀) are close to each other. The theoretically calculated bulk modulus (54 GPa) coincides well with the experimentally derived value (57 GPa) of the 3R phase.^[32] A notable feature is that the volume of the I4/mmm phase shows an abnormal change between 31 GPa and 32 GPa, which is assigned to the transition from the low pressure (LP) to high pressure (HP) phase [Fig. 3(d)]. The bulk modulus of the I4/mmm phase increases from 47.66 GPa to 135.07 GPa from LP to HP range. The HP I4/mmm phase is harder to be compressed than the other three layered phases (1T, 2H, and 3R phases) at ambient condition, which agrees well with the 3D structural characteristics of the HP I4/mmm phase [inset in Fig. 6(d)] with respect to the 2D layered structures of the 1T, 2H, and 3R phases (Fig. 1).

Fig. 3.

	Figure Option ViewDownloadNew Window
	Fig. 3. Unit cell volume variations with the increasing pressure for (a) 1T, (b) 2H, (c) 3R, and (d) I4/mmm phases. The solid lines are fitted by BM-EOS.

The energetic mechanism and kinetic mechanism are key factors to determine whether a phase can stably exist or undergo a transition under pressure.^[51–53] In order to explore the structure stability of NbS₂ under high pressure, we plot the enthalpy difference (ΔH) for these experimentally already-synthesized structures (1T, 2H, and 3R phases) relative to the theoretically predicted I4/mmm phase in Fig. 4 with vdW interaction of opt-B88. The previous studies on TiO₂ polymorphs indicate that the vdW interaction is an indispensable factor for describing the exact energy and enthalpy under pressure, so that the accurate phase transitions can be captured in accordance with the experiments.^[54] At ambient pressure condition (0 GPa), the enthalpy of the 2H phase is comparable to that of the 3R phase, which is lower than the 1T one. The theoretically calculated results are in line with the fact that the 1T phase is hard to be synthesized in bulk, only reported in thin film^[27] or monolayer forms.^[55] In contrast, the 2H and 3R phases are two common polymorphs been successfully synthesized in bulk forms.^[28,31]

Fig. 4.

	Figure Option ViewDownloadNew Window
	Fig. 4. Calculated enthalpy differences (ΔH) of the experimentally synthesized phases (1T, 2H, and 3R) relative to the theoretically predicted I4/mmm phase under hydrostatic pressure.

In addition, the enthalpy of the I4/mmm phase is far higher than those of the three experimental phases (1T, 2H, and 3R) at ambient pressure. These tendencies are consistent well with the previous reported theoretical results.^[37] As shown in Fig. 1, although the Nb atoms are six-coordinated with S atoms in all three experimental phases, the coordination polyhedra around the Nb atoms are different. The 1T phase contains hexagonal tightly packed S atoms and edge-sharing NbS₆ octahedra [Fig. 1(a)].^[27] Both the double-layer 2H and triple-layer 3R structures are consisted of trigonal prisms with different interlayer relative positions of the polyhedra [Figs. 1(b) and 1(c)], whereas the I4/mmm phase is a layered structure with face-sharing eight-coordinated body-center-cubic polyhedra [Fig. 1(d)]. When the pressure increases to 17 GPa, the six-coordinated 3R phase transforms to the 1T phase with the coordination polyhedra changing from triangular prisms to octahedra, but the coordination number keeps unchanged. Under further compression, a transition from 1T phase to I4/mmm phase takes place at 35 GPa, accompanying with a volume collapse of 8.1% [Fig. 5(a)]. Consistent with the previous experimental results that the 2H phase is stable up to 20 GPa,^[34] the 2H phase maintains thermodynamic stability up to 38 GPa under compression and then transforms to the I4/mmm phase, accompanying with a volume collapse of 7.8% [Fig. 5(b)]. The volume evolutions show discontinuous changes at the transition pressures, suggesting that the structural transitions are first order. It should be noted that the transition pressures in the present work differ from the previous theoretical results obtained by Liu et al.^[37] The differences mainly come from the choice of the exchange–correlation parameters (GGA vs. LDA) and the vdW corrections included during the geometrical optimization.

Fig. 5.

	Figure Option ViewDownloadNew Window
	Fig. 5. Variations of the volume as a function of pressure. (a) Structural phase transitions from 3R to 1T phase at 17 GPa, and from 1T to I4/mmm phase at 35 GPa. (b) Structural phase transition from 2H to I4/mmm phase at 38 GPa.

Another noteworthy point is that the slopes of the enthalpy increase significantly around the pressure range of 30–35 GPa (Fig. 4). In order to inspect the detail of the abnormal change of the enthalpy and the lattice volume variation of the I4/mmm phase, we decrease the pressure step to 1 GPa in the range of 30–35 GPa, and increase the pressure range up to 80 GPa [Fig. 3(d)]. Both the enthalpy and energy show abnormal changes between 31 GPa and 32 GPa [Figs. 6(a) and 6(b)]. The volume collapse causes a slight decrease of enthalpy and an abnormal increase of energy at 32 GPa. Similar to the cases of the 2H and 3R phases, the compression behavior of the I4/mmm phase shows strong anisotropy.^[32–34] The lattice constant c obviously decreases with increasing compression, whereas lattice constant a is almost unchanged. The compressibility is consistent with the layered structure stacking along the crystallographic c axis, arising from the weak interlayer vdW interactions. Along with the pressure increasing up to 32 GPa, the lattice constants show discontinuous changes, but the crystal symmetry and the space group are maintained. As shown in Fig. 6(c), the lattice constant a expands a bit and c shrinks a lot at 32 GPa; the lattice constants are a = b = 3.01 Å, c = 9.35 Å at 31 GPa and a = b = 3.17 Å, c = 7.99 Å at 32 GPa. The variations of enthalpy, energy, and lattice constants indicate that an isostructural transition occurs in the tetragonal I4/mmm phase.

Fig. 6.

	Figure Option ViewDownloadNew Window
	Fig. 6. Variations of the (a) enthalpy, (b) energy, (c) lattice constants, and (d) interlayer and intralayer Nb–S bond lengths as a function of pressure for the I4/mmm phase up to 80 GPa. The inset in (d) shows the structure of HP I4/mmm phase. The transition region is highlighted around 31–32 GPa.

In order to get further insights into the structural phase transition in the I4/mmm phase, we calculate the electronic structure of the I4/mmm phase around the pressure range of 31–32 GPa. The band structures and density of states (DOS) are displayed in Fig. 7. As shown in the band structures, both the conduction bands and valence bands overlap with each other near the Fermi level, displaying good metallic properties. The orbital overlap and atomic interaction are enhanced greatly under the compression. As shown in Fig. 7(a), the bands near the Fermi level show weak dispersion along the Z–R–X–Γ direction, leading to localized electron states with heavy effective mass. By contrast, the bandwidths are obviously broadened along the M–Γ–Z direction, leading to the delocalization of electrons. By further increasing the pressure, the out-of-plane dispersion of the bonding orbitals and antibonding orbitals increases significantly along the Z–A–M direction [Fig. 7(b)]. These changes of the band structures indicate an enhancement of bonding effect, and new bonds emerge in the crystal. As shown in the electronic structures, the Nb d states strongly hybridize with the S p states, and make major contributions to the bands around the Fermi level. Especially, a sharp peak locates close to the Fermi level in the DOS at 31 GPa, consistent with the almost degenerate band lying near the Fermi level along the Z–R direction [Fig. 7(a)]. The close proximity of this DOS peak to the Fermi level favors an electronic instability at 31 GPa.^[56] With further compression, a phase transition will take place to eliminate the instability. The sharp peak indeed disappears and the DOS at the Fermi level reduces at 32 GPa [Fig. 7(b)]. The electronic instability is close related to the interlayer vdW interactions. The interlayer Nb–S distance is 3.30 Å at 31 GPa [Fig. 6(d)], indicating that the interlayer vdW interaction is at the edge of existence (3.1-5.0 nm), causing an instability of electronic state. Further compression reduces the interlayer Nb–S distance to 2.75 Å at 32 GPa. Consequently, the interlayer distance is below the threshold of vdW interactions, which leads to that the interlayer vdW forces weaken and the electronic structure transforms to a new stable state.

Fig. 7.

	Figure Option ViewDownloadNew Window
	Fig. 7. Electronic band structure and density of states of I4/mmm-NbS2 at (a) 31 GPa and (b) 32 GPa. The colors of the band structure denote the contributions from Nb d (red) and S p (blue) states, and strong hybridizations between them.

The origin of the transition around 31 GPa can be further explored by the electron localization function (ELF).^[57] The ELF is a function ranging from 0 to 1, indicating the type and strength of the bonding. As shown in Fig. 8, the ELF is close to 0.9 around the S atoms, indicating that polar bonds exist. Electrons are transferred from Nb to S and molecular orbitals are formed. It is clearly seen that a layered structure exists at 31 GPa with eight-fold coordinated Nb [Fig. 8(a)]. However, the layered structure has become blurred at 32 GPa [Fig. 8(b)]. The vdW radii of Nb and S ions are 2.15 Å and 1.8 Å, whereas the covalent radii of Nb and S ions are 1.64 Å and 1.05 Å, respectively.^[58,59] The interlayer and intralayer Nb–S distances are 3.301 Å and 2.536 Å at 31 GPa, and 2.731 Å and 2.574 Å at 32 GPa, respectively [Fig. 6(c)]. Apparently, the interlayer Nb–S distance has reached the requirement of bonding at 32 GPa. New interlayer Nb–S bonds formation leads to the coordination number of Nb increasing from eight to ten, and the interlayer vdW interaction no longer plays the dominant role, but is replaced by a negligible weak interaction described by standard GGA.^[60] Accordingly, the I4/mmm structure transforms from a layered phase to 3D bulk phase due to the interlayer bonding effects, resulting in a higher compactness and lesser compressibility for the I4/mmm HP structure as shown in Fig. 6(d). Pressure-induced decrease or increase of the inter-atom distance usually causes variations of the coordination number. Similar phenomenon has been reported in Bi₂Te₃ layered structure, where the pressure breaks the layered structure with coordination number increasing from six to eight.^[61]

Fig. 8.

	Figure Option ViewDownloadNew Window
	Fig. 8. ELF projected on the (110) crystalline plane for the I4/mmm phase at (a) 31 GPa and (b) 32 GPa. The scale of ELF is ranging from 0 (blue) to 1 (red) and the interval is 0.1.

We further examine the structural stability of the I4/mmm phase under tensile strain in ab-plane. Accompanied with the in-plane stretching, the lattice constant c and the interlayer Nb–S distances are gradually compressed. As shown in Fig. 9, the lattice constants c, c/a ratio, and interlayer Nb–S distances show discontinuous changes with the tensile strain increasing from 8% to 9%. The energy evolution also shows different variation tendency around this strain range. With the tensile strain increasing up to 8%, the interlayer Nb–S distances gradually shrink from 4.337 Å to 3.696 Å. The lattice constant c compression accompanied with the ab-plane stretching renders the interlayer Nb–S distances on the verge of the vdW interactions. Further increasing the tensile strain, the interlayer Nb–S distances collapse and the interlayer vdW interactions disappear, concomitant with new covalent bond emerging and structure transforming from layered 2D to 3D bulk. The calculated electronic structures show similar character with the cases under high pressure, indicates that an instable electronic state becomes stable under compression. The interlayer bonding of Nb and S triggers the structural transition from 2D layered structure to 3D bulk phase. We can expect the structural transition tuned by the tensile strain to be realized in epitaxial films on appropriate substrates.

Fig. 9.

	Figure Option ViewDownloadNew Window
	Fig. 9. Variations of the (a) lattice parameters, (b) c/a ratio, (c) inter-atomic distances, and (d) energy under tensile strains in ab-plane of the I4/mmm phase.

In summary, the stabilities of NbS₂ under hydrostatic pressure have been investigated systematically based on DFT calculations with the inclusion of vdW corrections. We not only confirm the previous theoretically predicted structural transition, but also present electronic structures in better agreement with the available experimental results of the experimentally already-synthesized bulk phases. Furthermore, the structural transition from 2D layered structure to 3D bulk phase in the tetragonal I4/mmm phase is explored at 32 GPa, which can be achieved under tensile strain in ab-plane. The unusual structural transformation and electronic structure variation interplay with the change of interlayer bonding. The present work not only highlights the vdW interactions in forming the layered crystal structure, but also implies the critical role played by the vdW interactions in the structural transition of NbS₂. We provide new insights into the description of the structural stability of NbS₂ and complement the content of phase transition engineering.